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' A theoretical model is proposed to explain the roughness characteristics of an ice 

surface grown from a gravity and wind-driven supercooled water film flowing over 
Ch ■ 

K^ ' an inclined plane. The effects of the water supply rate, plane slope and air stream 

^ , 

^ ■ velocity on the spacing and height of ice surface roughness are investigated from 

a new type of morphological instability of the ice-water interface. The proposed 
macro-scale morphological instability under a supercooled water film is quite different 
from the micro-scale one which results in dendritic growth. It was found that ice 
, Ph| surface roughness spacing depends mainly on water layer thickness, and that surface 



roughness height is very sensitive to the convective heat transfer rate at the water- 
air interface. The present model takes into account the interaction between air and 
water flows through the boundary conditions at the water-air interface. This leads 
us to a major finding that tangential and normal shear stress disturbances due to 
airflow at the water-air interface play a crucial role not only on the convective heat 
transfer rate at the disturbed water-air interface but also on the height of the ice 
. surface roughness. This is confirmed by comparison of the amplification rate of the 

I ice-water interface disturbance predicted by the model with the roughness height 

observed experimentally. 
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I. INTRODUCTION 



A variety of surface features of growing crystal under a thin layer of moving fluid, 
which separates the developing solid from the surrounding air, are observed in natural 
phenomena.-'^. A first example is the ring-like ripples on the surface of icicles.^ A pattern 
similar to icicle ripples can be experimentally produced on the surface of a wooden round 
stick and that on a gutter on an inclined plane in a cold room, as shown in figures [H (a) 
and (b).-i^ These ripples appear clearly when water dripping from the top of the stick and 
gutter spreads effectively and covers the entire ice surface uniformly. The latent heat from 
the ice-water interface to the environment through the water film must be released during 
the freezing process. Consequently, a negative temperature gradient develops ahead of the 
growing ice beneath the water film and the water is a supercooled.- The spacing between 
ice ripples formed on the vertical stick and gutter was nearly 1-cm long, like natural icicle 
ripples. The wavelength of the ice ripples on the gutter decreases as the slope of the inclined 
plane increases. It increases only gradually as water supply rate increases, and the ripples 
move upwards very slightly with time.-i^ 

As the second example, figure [1] (c) is a schematic view of ice roughness formation on 
the parabolic leading edge of a NACA 0012 airfoil under glaze icing conditions (i.e. air 
temperature close to freezing and high liquid water content (LWC)), observed by Shin.- 
LWC is the mass of water contained in a unit volume of air. In glaze icing, the portion of 
the impinging water droplets that cannot be frozen runs off the surface due to gravity or 
wind drag. The latent heat released in the freezing process must be transferred from the 
ice-water interface through the unfrozen water film to the air.-i2 Shin defines roughness as 
surface irregularity growing on the top of macro-ice shape with horns and feathers. Smooth 
to rough zones in figure [1] (c) is defined as the region where surface condition changes 
from a smooth to a rough one. Roughness size was measured for various airspeeds, air 
temperatures and LWCs. Roughness height increases with increasing air temperature and 
LWC, whereas airspeed has little effect on roughness height. Roughness spacing is of the 
order of a millimeter, decreases with increasing airspeed, and increases with increasing air 
temperature and LWC. The boundary between smooth and rough zones moves upstream 
towards stagnation region with time, as shown in figure [T] (c). 

For the third example, figured] (d) is a schematic view of an initial aufeis (also referred to 
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as icings) formation in frigid air when a shallow sheet of water, introduced at the upstream 
end of the wind tunnel, flows or trickles over a sloped frigid surface.— Both gravity and wind 
drag drive the spreading of shallow flows of freezing water. The initial aufeis morphologies, 
characterized by wavelike or terraced forms, are shown in figure [T] (d). Their roughness 
spacing and height was found to vary with slope and wind speed. As the slope increases, 
their roughness spacing decreases and roughness height increases. Moreover, their roughness 
spacing and height decrease as wind speed increases.— 

Finally, travertine terracing is among the most spectacular geological phenomena on 
earth, not only in limestone caves and around hot springs, but also in streams and rivers in 
limestone terrain. The interactions between hydrodynamics, water chemistry, calcium car- 
bonate precipitation and carbon dioxide degassing constitute a complex pattern formation of 
travertine terracing.— The relationships between slope, discharge, terrace wavelength, depth 
and height are discussed by Pentecost™ Inter-dam distances increase with large discharge, 
where dam is defined as terraces that are filled with water, forming pools and lakes. Pools 
are shorter on steep slopes, instead height is larger. Interestingly, variations of the inter-dam 



distances and ice ripp 
and Fig. 8 (a) in Ref. 
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e wavelength with slope show the same trend (see Fig. 5 in Ref. 
sl), regardless of different crystals. 
We can see common features for the surface roughness characteristics such as roughness 
spacing and height among phenomena mentioned above. For the first example, a theo- 
retical model of the origin of ripples on icicles was proposed,—"— and the results were in 
good agreements with the experimental results.-"^ For the second example, in order to ex- 
plain experimental results on glaze ice roughness diameters, accreted on NACA 0012 airfoil 
leading edges reported by Shirt^, Tsao and Rothmayer developed a high Reynolds number 
triple-deck theoryi^ to describe the interaction between the air boundary layer, water film 
and glaze ice sheet.— A novel broad-band ice instability mode was found in regimes with 
simultaneous air and wall cooling, but there was no well-defined maximum amplification 
rate wavenumber (or equivalently, wavelength). To overcome this issue, the Gibbs-Thomson 
effect was introduced to stabilize the smallest scale icing disturbances. However, the length 
scale predicted by their theory was much smaller than the roughness spacing of the order 
of millimeters observed in the experiments of Shin.^^ For the third and final examples, the 
quantitative morphology and size distribution of aufeis and travertine terraces as a function 
of parameters such as slope, water fiux and airspeed have not yet been studied in detail. 



^ V ^ Impin&ing droplets 




height 



FIG. 1. (a) Ice ripples on (a) a stick and (b) a plane.- (c) Schematic of ice roughness formation 
on the parabolic leading-edge of a NACA 0012 airfoil in glaze icing conditions.^ (d) Schematic of 
an initial aufeis formation on a sloped surface.— 
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Therefore, it is necessary to construct a comprehensive model to elucidate common features 
for the roughness characteristics described above. 

In this paper, in [TTl we propose a theoretical model to explain the effects of water supply 
rate, plane slope and wind speed on the roughness spacing and height of an initial aufeis 
(icings), from the morphological instability of growing crystal under gravity and air shear 
stress driven water film. In a previous ice growth model, water is driven by gravity alone.- 
In another model, water is driven by air shear stress alone,— but the model is valid on 
a horizontal surface. However, by combining the two driving forces, the resulting model 
becomes more complex than the previous ones because air and water flows and temperature 
fields are highly coupled with the water film thickness. Therefore, a numerical method 
is proposed to solve the governing equations for an air-water-ice multi-phase system. In 
llllt the experimental results concerning the roughness characteristics of the initial aufeis 
observed by Streitz and Ettemar^ will be explained theoretically. It will also be shown 
that the growth conditions of the ice-water interface disturbances are strongly affected by 
variable air stresses exerted on the water-air interface by the airflow. Crucial evidence of the 
importance of such air shear stress disturbances will be shown by comparing theoretically 
calculated amplification rates of the ice-water interface disturbance with the ice surface 
roughness heights observed by Streitz and Ettema. Concluding remarks are made in IIVI 

II. MODEL 

In the experiments of Streitz and Ettema^° for an initial formation of aufeis on an inclined 
plane shown in figure [1] (d), a wind tunnel set up in a refrigerated laboratory was designed 
to tilt only downwards, allowing the water to flow as a thin sheet driven by gravity and 
wind drag. The water originates from a row of holes located at the top of the plane. The 
experiments was conducted with a water flow rate of about 1692 [(ml/h)/cm] at an initial 
water temperature of 0.1 °C, in air at a temperature of -5 °C. The initial formation of aufeis 
was defined as the initial layer of aufeis formed when shallow water first flows as a laminar 
sheet over a frigid surface and freezes onto it. The initial morphologies appeared essentially 
wavelike or terraced, and their spacing and height indicated in figure [T] (c) were measured 
for plane slopes up to 15° and wind speeds up to 48 km/h. 

The current model configuration and coordinate system is shown in figure [2l which is 
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FIG. 2. Diagram of physical model and coordinate system (vertical height is not to scale). The solid 
curves represent the locations of the disturbed ice- water and water-air interfaces. The dashed curves 
represent the locations of the undisturbed water-air interface and air boundary layer thickness. 

based on the laboratory experiments of Streitz and Ettema mentioned above, x is the 
position along the inclined plane measured from the location of the water source, and y is 
the position measured from a fiat ice- water interface. 6 is the angle of the inclined plane with 
respect to the horizontal. Air temperature is at Too, which is lower than the temperature 
Tsi at a fiat ice-water interface. We assume a steady laminar airflow parallel to the x axis 
with free stream velocity Uoo- By releasing the latent heat to the air through the water-air 
interface at temperature T;^, ice grows from the portion of the supercooled water film driven 
by gravity and wind drag force exerted by the airflow. 6 and Jiq are the air boundary layer 
and water film thickness, respectively, and uia is the water surface velocity. As shown in 
figure [2], the water-air interface and flow and temperature in the air are disturbed due to 
a change in ice shape. As a result, the air shear stress exerted on the water-air interface 
and the heat transfer rate from the water-air interface to the air are variable. These in turn 
affect the flow and temperature distributions in the water layer. In response to this change, 
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further ice-water interface development is determined. In this sense, the air boundary layer, 
the unfrozen water and the ice substrate become a complex air-water-ice multi-phase system. 

The following assumptions are used in the present model: (1) Water is discharged from 
only the top of the plane and there is no airborne water droplets impinging on the plane 
surface. (2) The free stream velocity Uoo is constant in space. (3) Density remains constant 
through the phase change. (4) Due to the long time scale of the ice-water interface motion, a 
quasi- stationary approximation is used for the disturbed fields, and unsteadiness only enters 
through the Stephan condition. (5) Even in the presence of the undisturbed temperature 
gradient in ice, the morphological instability occurs when the ice thickness exceeds a critical 
thickness. As far as the ice thickness is large, it is a good approximation to neglect heat 
conduction into a substrate beneath the ice sheet.— (6) The presence of waves on the water 
film is ignored because the waves did not interact with the forming ice in any observable 
manner in the experiments, except for enhancing the spreading of the water over the aufeis 
surface.— (7) Freshwater icing sponginess containing non-negligible amount of liquid water 
was observed in aufeis.— The spongy ice formation, in which a portion of the surface liquid 
is incorporated into the ice matrix, is not considered. 



A. Governing equations 

The following basic equations and boundary conditions governing the air-water-ice multi- 
phase system are based on previous papers^"^ and are reviewed here to ensure a relatively 
self-contained treatment. The velocity components in the x and y directions in the air, Ua 
and Va, are governed by 

dUa , dUa , dUa 1 dpa , f d'^Ua , d'^Ua\ ... 
+ Ua^— + Va^- = ^ + l^a + , (1) 



ot ox oy Pa ox \ Ox'^ oy 



dt dx dy Pa dy " \ dx'^ dy"^ J ' 

dua dva „ ,„s 
ox oy 

where pa is the air pressure, pa = 1.3 kg/m^, the density of air, and z/q = 1.3 x 10"^ m^/s, 
the kinematic viscosity of air. The velocity components in the x and y directions in the 



water layer, ui and vi, are governed by 



^ j^u — +v — -- — — + u + ] + sine 

dt dx dy pi dx \ dx"^ dy'^ ' 



(4) 



dvi ^ dvi ^ dvi 
dt dx dy 



1 dpi ( d'^vi d'^vi\ _ 
+ + ~ gcosO, 

pi oy \ ox^ oy"^ 



dui dvi 

— - H ^ = 0, 

dx dy 



(5) 



(6) 



where vi = 1.8 x 10~^ m^/s and pi = 1.0 x 10'^ kg/m^ are the kinematic viscosity and density 
of water, respectively, pi is the water pressure, g is the gravitational acceleration, and 6 
is the plane angle. The continuity equations (|3]) and ([6]) can be satisfied by introducing 
the stream functions ipa and ipi such that Ua = dipa/dy, Va = —dipaldx^ ui = dipi/dy, and 
vi = —dipi/dx. 

Both velocity components ui and vi at a disturbed ice- water interface, y = ((t,x), must 
satisfy the no-slip condition: 



Ul\y = ^ = 0, Vl\y = (^=0. 



(7) 



Since there is no impingement of supercooled water droplets on the water film, the kinematic 
condition at a disturbed water-air interface, y = ^(t,x), is 

de , d^ 



The continuity of velocities of water film fiow and airfiow at the water-air interface is 



The continuity of tangential and normal stresses at the water-air interface is 



(8) 



(9) 



/ dui 




dvi 




[ dUa 




^ dVa 




V 9y 


H 


dx 






y=^ 


dx 


y=i) 



-Pa\y=ii + 



dVg 

dy 





-Pi\ 







dvi 
dy 



y=i 



-7 



dx'^ 



1 + 



dx 



-3/2 



(10) 



where pi = piUi = 1.8 x 10 Ns/m^ and pa = Pa^a = 1-69 x 10 ^ Ns/m^ are the viscosities 
of water and air, respectively, and 7 = 7.6 x 10~^ N/m is the water-air surface tension. 
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The equations for the temperatures in the air T^, water Ti and ice Tg are 

m OTa dTa fd^Ta d^Ta\ 

+ Ua^T- + Va^T- = Ka — + — , (12) 



01 OX ay \ ox"^ oy 

— +U — + V—-K + (13) 

dt dx dy ^ \ dx'^ dy"^ ) ' 

^ + (14) 



dt \ dx'^ dy"^ y 

where Ka = 1.87 x 10~^ m^/s, ki = 1.33 x lO^'' m^/s and Ks = 1.15 x 10~^ m^/s are the 

thermal diffusivities of air, water and ice, respectively. 

The continuity condition of temperature at the ice-water interface is 

Tl\y=t; = Ts\y=(; = Ti, (15) 

in which the interfacial temperature Tj is an unknown to be determined. The conventional 
Stefan problem cannot describe the pattern formation observed in nature.— Likewise, all 
linear stability analyses based on the assumption that the temperature at the disturbed 
ice-water interface remains at the equilibrium freezing temperature, = Ts\y=(^ = Tgi 

{Tgi =0 °C for pure water), showed that the ice-water interface disturbance becomes unstable 
for all wavenumbers, and that there is no dominant amplification rate to select a preferred 
wavelength.-i^'^ The Stephan condition is 



y=c. oy 



(16) 

s/=C 



dt J dy 

where L = 3.3 x 10^ J/m^ is the latent heat per unit volume, V is the undisturbed ice 
growth rate, and Kg = 2.22 J/(mKs) and Ki = 0.56 J/(mKs) are thermal conductivities 
of ice and water, respectively. 

The continuity condition of temperature at the water-air interface is 

Tl\y=^ = Ta\y=^ = Tla, (17) 

where Tia is a temperature at the water-air interface and will be determined later. The 
continuity of heat flux at the water-air interface is 

3r„ 



dy 



y=<i dy 



(18) 



where Ka = 0.024 J/(mKs) is the thermal conductivity of air. Far away from the air 
boundary layer, the velocities and temperature asymptote to their far-field values: 

^a|j/=oo ^oo) ^a|j/=oo 0; -^a|j/=oo -^oo. (1"^) 
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B. Linear stability analysis 



Since most of the derivations of the stability analysis follow the same procedure given in 



Refs. 



and 



2(1, except for the modified undisturbed velocity profile in the water film and 
the boundary conditions of ([25]), ([26D, ([36]), §7^, (gO]), dH]) herein, full details will not be 
given. 



The field variables in § 111 Al are assumed to be decomposed into undisturbed and disturbed 
parts, as follows: 







1 o\ 










ho 










^a 




Uoofaijfjik 


Ipl 




'4'1 




Ulafl{y*)Ck 


Pa 




Pa 


+ 


{paUl^/So)ga{vKk 


Pi 




Pi 




{piu1a/ho)gi{y*)Ck 


Ta 




fa 






Tl 




fl 






[tJ 








\ Hs{y*)GiCk ) 



exp[(Tt + ikx\. 



(20) 



A simple normal-mode analysis is applied to the ice-water interface disturbance C, and the 
corresponding fields variables ,%lj'^,%lj[,p'^,p'i,T'^,Tl,Tl), which are the disturbed part in 
( l20l) . Here ho is the undisturbed water film thickness, uia is the surface velocity of the water 
film driven by gravity and air shear stress, 5o = {2h'ax/uooY^'^ is a scaled measure in the 
air,— Moo is the free stream velocity, x is the distance from the leading edge where water is 
supplied, r] = (y- ho)/5o, y* = y/ho, Ga = -dfa/ dy\y^i^, Gi = -dfi/dy\y=o, and (k and 
are the amplitudes of the ice- water interface and water-air interface, respectively. ga, Ha) 
and {fl, gi, Hi, Hg) are dimensionless functions with respect to r] and y*, respectively, k is the 
wavenumber and a = a^^'^ +ia^^\ a^^^ and Vp = ~a^^^ /k are the amplification rate and phase 
velocity of the disturbance, respectively. Since the undisturbed part of heat conduction in 
the ice is assumed to be zero, Tg = T^i is used throughout this paper. 
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1 . Governing equations for undisturbed and disturbed parts of flow and 
temperature in the air 



Substituting ipa = i^a + i^'a = UooSoFa{r]) + u^fa{rj)CkewWt + ikx] andT„ = T^ + T^ = Ta + 
Ha{r])Ga^kGxp[at + ikx] into the complete equations ([I]), ([2]) and ( fT2l) . a set of dimensionless 
differential equations for the undisturbed part Fa, = (Ta — Too)/ (T/a — Too) and for the 
disturbed part fa, Ha are obtained: 

d'Fa ^ d^F„ 



-Fa- 



dr] 



a 
2 ' 



dri^ 



d'fa 



dif 



+ <^ 2kl^ - (2 - iK^Rta) 



d'^T dT 

' a^^ a' 



drj 

dFa\ d^fa 



(21) 
(22) 



di] J dr]'^ 



+ Hl{Fa + 27]- 



dFa\ d-'FAdf, 



di] 



dr]"^ J dr] 



dF„ d^F„ 



+ 



drj dr]^ 



fa. 



(23) 



d'^{Ga*Ha) 



-Pr F 



d{Ga*Ha, 



drj 



+ \kl, + Pra{-1 + ika^Rea)^^ [GaMa] 



-ika^iPraRCa 



dTa. 
drj 



J a, 



(24) 



where Rca = u^So/va and Pr^ = Va/i^a are the Reynolds number and the Prandtl number 
of air, respectively, ka* = k6o is the dimensionless wavenumber normalized by the length 6q, 

and Ga* = -dfa*/dT]\r,=o. 



The undisturbed part of (19]), Ua 
Too, the disturbed part of ([9]), u'^ 
disturbed part of f lT7|) and T^ 
fa and i^a, respectively, 



y=oo = dipa/dyly 
y=oo = d^j'a/dyly 

y=OD 



Tin and T 



= -dip'a/dx\y 

yield the following boundary conditions for Fa, Ta*, 



= v'\ 



a\y=OD 

0, the 



dF„ 



drj 



_ Uia_ 
77=0 Uoo ' 



a\ri=0 



0, 



c?T„ 



- a* I rj=0 



dfa 

drj 



r;=0 



drj 

^ \ =0 

-a*\ri=oo 

So Uia dfi 



rj=oo 



(25) 



p I Ula J, I 

Ja r7=0 // ■(/*=1' 

7/ 



r?=o /lo Moo dy* 

dfa 



drj 



0, 



/ fl\y,^ 
y,=ll 

fa\ 



r)=oo 



I r]=oo 



H, 



a\ri=0 



H. 



a I r]=oo 



(26) 
(27) 
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2. Governing equations for undisturbed and disturbed parts of flow and 
temperature in the water film 



The undisturbed part of ffTSl) yieldsi^i*^ 

K h 

Substituting ( l28l) into the undisturbed part of ( |T6i) . the undisturbed ice growth rate is 
obtained:-^'20 

K T 

V = ■ (29) 

L{6o/Ga*) 

The length 60/ G a* in ( |28|) and ( 129|) is regarded as the air boundary layer thickness 6 in 
figure [2j 

Substituting ^Ji = ipi + ipl = uiahoFi{y^) + uiafi{y*)(kewWt + ikx] and Ti = Ti + TI = T; + 
Hi{y:^)Gi(kG:^v[o't + ikx] into the complete equations (jl]), ([5]) and ( fT3|) . a set of dimensionless 
differential equations for the undisturbed part ui^ = ui/uia = dFi/dy^, Ti^ = (Ti — Tsi)/{Tsi — 
Tia) and for the disturbed part //, Hi are obtained: 



dyl viuia 
d'^fu 



(30) 



^ = (2fcf* + ikuReiUu) ^ - j/cf^ + ifc^^-Re^ + ^p^^ | /z, (32) 

d ^^1*^1) ^ ^/.2^ _^ ikuPeiUu) [GuHi] - ik^Pei^^fi, (33) 

where Rei = uiJiq/i'i and Pe; = uiJiq/ki are the Reynolds number and Peclet number of 
water, respectively, fc/* = kfiQ is the dimensionless wavenumber normalized by the length Jiq, 
and Gu = -dfu/dy^\y^=o. When deriving ho/iTsi - Tia)d{Tsi - Tia)/dx = dho/dx - 
hQ/{2x) is used, which is obtained by differentiating fl28|) with respect to x. Equations (|30|) - 
(13 3 p are finally obtained by neglecting the term with dJiQ/dx <^ 1 and Hq/x <^ 1. This is in 
agreement with the more usual lubrication approach .-"^^i^ Using the boundary conditions 
ui\y=o = 0, nidui/dyly^h^^ = Hadua/dyly^h^, 7]|y=o = T^i and fi\y^-f^^ = Tia, the solutions of 
the dimensionless undisturbed velocity and temperature profiles in the water film are given 
by 



ghl sin 6* 2 f 9^1 sin 9 /iaWoo^o d'^^a 



ui* = — 7; y* + ^ 



^l^lUia * V '^I'^la l^mJo dlf 
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7?=0 



y*, Tu = -y*. (34) 



Linearization of the disturbed part of ([7]) at y = 0, as well as (fTOj) . (fTT]) . ( !T7|) and (fT8|) at 
y = ho yield the boundary conditions for /; and if/: 

dy* 



0, 



dyl 



d'^uu 



dyl 



fl\y,=0 — 
Sal// 



/ls/.=i - 

- (Sfcf^ + ikuRei) 
y.=i ay* 



(35) 
(36) 



+ikuRei 



duu 



dy* 



dyl 



= 0, 



5o 



J7=0 



Here, in (l36l) and lia in (1371) are defined by 



^0 



Pl \Ul 



Pl Ula 
2 



X 







drf 


ri=0 { 



hp 1 + ikg^Reg 
5o 1 + {ka^RtaY 

- Ukl^ + {-I + ika.Rea 



0, 

+ ^a*/a|)?=0 



(37) 
(38) 
(39) 

(40) 



dFa 



drj 



1 ^ 

r]=Q J (i?7 



77=0 



-\-ika*Rc, 



dif 



r;=0 



/a I r;=0 



(41) 



where Fr = uia/{ghoY^'^ is the Froude number, and We = 'j/ipiufjio) is the Weber number. 
It should be noted that the disturbed part of the water flow is affected by the airflow through 
the terms in (1361) and Ua in fl37l) . which are hereafter referred to as the tangential and 
normal air shear stress disturbances, respectively. 

Using fl34l) and u'l = dip'Jdy , the volumetric water flow rate per width is given by 



Qlk 



[ui + Ui)dy=- 



C 



3ui 



2niSo dr]'^ 



ri=0 



+uiahoi^k + fi\y,=iCk) exp[at + ikx]. (42) 
From the disturbed part of ([8]), the relation between the amplitude of the water-air interface. 



C,k, and that of the ice-water interface, is obtained: = —fi\y,=iCk- Therefore, 
be written as 



Q/L 



qsinO- o 
+ 



pa'^oo d Fa 



2/i;(5o dr]"^ 



r)=0 



hi 



can 



(43) 
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Applying the definition ui\y^ji^ 
is determined from 

Ula 



Ula (or equivalently, 1) to ([31]), tlie value of uia 



5f sin 61 -2 flaUood^Fa 



-hi + 



r?=0 



ho- 



rn 



2^1 ' fildo drf 

Linearizing the temperature at the ice- water interface in f|T5|) . Tj can be written as 
Tj = Tsi + ATs/, where Tgi is the temperature at an undisturbed ice- water interface and 
ATs; is a deviation from it when the ice-water interface is disturbed. Substituting T/ = 
-f^i(2/*)^iCfcexp[(Tt + ikx] and T'^ = Hs{y^)Gi(k^xp[at + ikx] into the disturbed part of 
( ITSjl and ( fT6l) . the dimensionless temperature deviation at the ice- water interface, ATgu = 
lm[ATsi/{Tsi — Tia)], the dimensionless amplification rate, ai^^ = a^^''^ /{V /Jiq), the dimen- 
sionless phase velocity, fp* = —a^^y^kV) are determined as follows: 



ATsu = Sb(t^) |(if/'"^|j^^=o - 1) sm[ki^{x^ - VpJ^)] + if/*^|y,=o cos[fc;*(a;* - VpJ^)]^ , (45) 



(r) 



(r) 



+ K^kUHt^\y,=o-l) 



y*=o 




dH 



dy* 



+ KfkuHP\y,=o 



y,=o 



(46) 
(47) 



where Im denotes the imaginary part of its argument, Hi"^^ and ij/*"* are the real and imag- 
inary parts of Hi, and Ki = K^/Ki = 3.96 is the ratio of the thermal conductivity of ice 
to that of water, x* = x/ho, t* = (y/ho)t, (5b(t*) = exp{ai^h^,)6b and 6b = Cfc/^o- It should 
be noted that ATsu — in the limit ki^ — 0, and AT,;* for a finite ki^ varies because the 
disturbed temperature distribution in the water layer. Hi, is affected by both air and water 



flows (see Ref. 



201 for more details). 



3. Numerical procedure 



Since the airflow was not considered in a previous paper ,^ Hq is determined by the gravity- 
driven part in (H^ . {Q/lw)g = (7^0 sin 6*/ (3z/;). This yields 

ft„= J^i'ji] (48) 



5'sin6' \/, 



and then 



and Ula = gh^ sin 6 / {2ui) from 
from (!34l) . Hence, the values dui^ / dy^\y^=o = 2, dui^/dy^ 



-y1 + 2y=K is the half-parabolic form 



and d'^uu/dyl\y^=i 
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are obtained. Noting that = in (l36l) and Hq = in (137|) in the absence of airflow, 
and that the term ki^Rei{cos9/ Fr"^ + WekfJ in fl371) is equivalent to the parameter a of 
(22) in Ref. is, the boundary conditions (1531) (15B]) and (^71) reduce to (23) in the previous 
paper." a was the parameter characterizing the effect of gravity and surface tension on the 
water-air surface. Furthermore, in the absence of airflow herein can be written as 
(PHa/drf = kl^Ha- Its solution is Ha = e~^°-*^ with the boundary conditions (1271) . which 
yields Jiq/ 6o{—dHa/ dri\rj=o) = ki^. Hence, the boundary condition fl39|) reduces to (33) in the 
previous paper.- 

In another previous paper,— since water flows on a horizontal ice surface, Iiq is deter- 
mined by the shear-driven part in fH51) . {Q/lw)s = fJ'aUood'^ F'a/ dri'^\r,=ohl/ {2fii6o) . Jiq can be 
expressed as 



and Uia = f^aUood'^Fa/dri'^\ri=oho/ ifiiSo) from (jH]), and then ui^ = y^. is the linear form from 
Hence, the values duu/ dy^\y^=Q = 1, dui^/dy^\y^=i = 1 and d'^ui^/ dyl\y^=i = are 
obtained. 

As Moo increases, the ratio of uia to Uoo approaches 0, as shown in figured (a) for various 
Q/lw Then the first equation of (!25|) and the first and second equations of (126|) can be 
approximated as dFa/dri\ri=o = 0, dfa/drj\r^=Q = —d'^Fa/drf[r^=Q and /a|r,=o = 0, using the 
fact that the viscosity ratio of air to water is very small, fia/fJ'i ^ 1- These conditions are 
equivalent to Ua\y=ho ~ ""£11^=^0 ~ ^ '^a\y=ho = O5 respectively, which means that the 
air effectively sees the water as a rigid body. Accordingly, the boundary conditions (l35l) . (l36ll 
and (Ell) with and dH]) reduce to boundary conditions (39), (40) and (41) with (42) and 
(43) in the previous paper.— Furthermore, as Mqo increases, the values of d'^Fa/ drf[r^=Q and 
Ga* = —dTa^l drj\^=Q converge to 0.47 and 0.41, respectively, for various Q/lw, as shown in 
figures [3] (b) and (c). Then the profiles Ua* and T^* are independent of the parameters Q/lwi 
9, Uoo and x and become similarity solutions for large Uoo-— When uia/uoo <^ 1, the solutions 
in the air are determined independently of the solutions in the water film. Hence, once Fa 
is obtained, ho is determined from (H9l) . On the other hand, in order to solve the governing 
equations for the water flow, the solutions of the airflow are necessary, as indicated in the 
terms and H^. 




(49) 
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FIG. 3. Variation of (a) uia/uoo, (b) d?Fa/dr]'^\rj=Q and (c) Ga* = —dTa*/dr]\r^=o with Uoo for 
Q/L = 1692(.), 1000 (■), 500 (♦) and 160/3 (a) [(ml/h)/cm]. 

By combining the two driving forces, the situation becomes more complex, as explained 
below. The air- water-ice multi- phase system considered here consists of ( 12T|) - (!24|) with 
boundary conditions f l2^ -f l27|) in the air, and (ES]) and fl55]) with and the boundary 
conditions in the water film, as well as a cubic equation fH5]) for Jiq. For a given 

Q/lw, 6, Uoo and x, the value of ho is numerically determined from fH3|) . However, the value 
dPFal drf\r^=Q is needed. When u^o is less than about 1 m/s, the term uia/uoo in f l25|) and fl26|) 
cannot be neglected. Furthermore, the first equation in (126|1 includes dfi/dy^\y^=i/ fi\y^=i, 
whereas Sa in ( l36i) and Ila in ( 1371) need the solutions Fa and /a as indicated in (HOj) and 
(HTl) . Since the solutions in the air and in the water film and ho are coupled, it is impossible 
to solve the current system simultaneously. 

To overcome this difficulty, the system is solved with the following iterative method. First, 
the temporal values of ho and dfi/dy^\y^=i/fi\y^=i are set. Then, the system is solved, except 
for (H3l) and the first equation in (!26l) . From these solutions and the excluded equations, new 
values of ho and dfi/ dy^\y^:=i/ fi\y^=i are obtained. Then, the system is solved again. After 
some iterations, ho and dfi/dy^\y^=i/ fi\y^=i settle to constant values, and the solutions of 
whole system is finally obtained. Substituting solution Hi into (H6i) and (H7|) and replacing 
k* by iho/So)ka*, cri''^ and Vp^ are presented with respect to ka*. 

The variation of the wavelength A of ice ripples formed on an inclined plane for various 
slope angles 6 of Fig. 8 in the previous paper- was obtained in the absence of airflow. It was 
confirmed that the current system reproduced the same variation of A with 6 at the limit 
of low wind speed, for example, Uoo = 0.01 (m/s), and for the same water supply rate of 
Q/lw = 160/3 [(ml/h)/cm] as used in the previous paper.- The case of u^o = is excluded 

16 



in the current system because Sq = {2i'aX /uaoY^'^ diverges. 
III. RESULTS AND DISCUSSION 

A. Variation of thickness Jiq and surface velocity uia with 9, Uoo and Q/lw 

Figures m (a), (b) and (c) show the variation of undisturbed water film thickness Tiq with 
9, Uoo and Q/lw, respectively, ho decreases with increasing 9 and Uoo, whereas it increases 
with Q/lw Figures m (d), (e) and (f) show the variation of the surface velocity of the water 
film, uia, with 9, Uoo and Q/lw, respectively, uia increases with increasing 6*, Uoo and Q/lw- 
Figure H] (a) shows that the differences between Jiq for different values of become smaller 
with 6', and that the asymptotic form of Jiq can be expressed by (jUj). This means that as 
9 increases, the water flow rate is dominated by the gravity- driven part in ( l43l) . {Q/lw)g, as 
shown by the solid curves in flgure H] (g) . 

On the other hand, in flgure H] (b), the differences between ho for different values of 9 
become smaller with Uoo, and the water flow rate is dominated by the shear-driven part 
in fH3|) . {Q/lw)s, as shown by the dashed curves in flgure S] (h). The asymptotic form of 
ho can be expressed by (149|) . From {Q/lw)g = {Q/lw)s, the shear-driven flow changes to a 
gravity-driven one at 



as 9 increases in flgure H] (g). On the other hand, the gravity-driven flow changes to a 
shear-driven one at 



as Moo increases, as shown in flgure H] (h). The points 9c and Mqoc move to the right in 
flgures m (g) and (h) as Uoo and 9 increase. 

B. Variation of wavelength A and phase velocity Vp^^ with 9, Uoo and Q/lw 

Figures|5](a) shows the variation of the dimensionless ampliflcation rates at'' = cr*^'''' / {V /ho) 
against the dimensionless wavenumber ka* = k6o for = 0°, 20° at u^o = 16 km/h and 
Q/lw = 1692 [(ml/h)/cm]. On the other hand, flgure[5](b) shows the variation of ai^^ against 
ka* for Uoo = 5, 30 m/s at 6' = 8° and Q/lw = 1692 [(ml/h)/cm]. An ice pattern with a wave 




(50) 




(51) 
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FIG. 4. For Q/l^ = 1692 [(ml/h)/cm]=4.7 x 10 ^(m^/s) and x = 0.1 m, variation of ho and 
with (a), (d) 9 for Uoo = 16,32,48 (km/h) (solid curves) and Uoo = 0.01 (m/s) (dashed curve); (b), 
(e) Uoo for 6* = 1°, 8°, 15°; (c), (f) Q/l^ for 6* = 8° and Uoo = 16 (km/h). Variation of gravity- 
driven water flow rate {Q/lw)g and shear-driven water flow rate {Q/lw)s with (g) 6 for Uoo = 16, 
48 (km/h) and (h) u^o for = 1°, 8°. 

number at which the amplification rate acquires a maximum is expected to be observed. 
For example, for 6 = 20° in figure O (a), ai^^ acquires a maximum value of aimax = 216.9 
at = 0.41. Since the wave number k is normalized by Sq, the corresponding wavelength 
of the ice pattern is 1.19 cm from A = 2n5o/ka*. Here, 5o = {2i'ax/uooY^'^ = 7.65 x 10~^ m 
estimated from x = 0.1 m and Mqo = 16 km/h is used. 

Using this method, the wavelength A was calculated for various 9 and Uoo- The variation 
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FIG. 5. (a) and (b) represent the variation of dimensionless amplification rate = a^^^ /(V /Tiq) 
against dimensionless wave number ka^, = k6o. Solid and dashed curves in (c) and (d) represent 
the variation of theoretically obtained wavelength A with 9 and Uoo- The symbol □ in (c) and 
(d) represent the experimental results of Ref. loj for the variation of roughness spacing with 6 at 
Uoo = 16 km/h and with Uoo at 6 = 8°, respectively, (e) and (f) represent the variation of A with 9 
and Uoo, respectively, for Q/lw = 500, 1000, 1692 [(ml/h)/cm]. (g) and (h) represent the variation 
of the dimensionless phase velocity = Vp/V with 6 and Uoo, respectively. The theoretical curves 
in (a)-(d), (g) and (h) are obtained for the same water supply rate of Q/l^ = 1692 [(ml/h)/cm] as 



used in the experiments of Ref. 



IC 



of A with 9 at = 16, 48 km/h and Uoo = 0.01 m/s is presented in figure |5] (c), whereas 
figure (d) shows the variation of A with Uoo at 6' = 1°, 8° and 15°. Here the theoretical 
results are shown by the solid and dashed curves. The symbol □ in Figs. M (c) and (d) 
represents the measured roughness spacing in the experimental for Uoo = 16 km/h and 
^ = 8°, respectively. Both theoretical and experimental results at Uoo = 16 km/h in Figs. 
|5] (c) show that A decreases with increasing 6. The variation of A for Uoo = 0.01 m/s also 
shows the same trend as that for Uoo = 16 and 48 km/h, despite large difference of Uoo- 
The experimental results (□) for = 8° in figure E] (d) shows that A rapidly decreases with 
increasing Uoo- On the other hand, theoretically obtained A for the same ^ = 8° gradually 
decreases for small values of Uoo, and increases very slightly with Uoo- For much larger values 
of Moo, which is not shown in figure |5] (d), A decreases again. The variation of A for 6* = 15° 
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also shows the same trend as that for 6 = 8°. X at 6 = 1° rapidly decreases for small values 
of Uqq and gradually decreases for larger Moo- Figures [5] (e) and (f) show the variation of A 
with 9 at Moo = 16 km/h and with Uoo at 6* = 8°, respectively, for Q/lw = 500, 1000, 1692 
[(ml/h)/cm]. These figures show that A increases with Q/lw It was found by comparing 
figures m (a) andO (c), as well as figures S] (b) andO (d) that the variations of A with 6 and 
Moo show almost the same trends as for Iiq. This indicates that Jiq is the most important 
parameter to determine ice roughness spacing. Since the experimental data (□) for 9 = %° 
in figure [5] (d) are scarce, it is difficult to conclude that there is large disagreement between 
the experimental and theoretical results. In addition, based on theoretical considerations, 
it seems impossible that the experimental results (□) shown in figure (d) decrease rapidly 
with Moo- In order for that to be true, Jiq aX 6 = 8° in figure S] (b) must decrease rapidly 
with Moo, like Tiq a,t 9 = 1°. 



Finally, figures [5] (g) and (h) show the variation of the dimensionless phase velocity 
Vp* = Vp/V with 9 at Moo = 16 km/h and with at = 8°. The magnitude of Mp* was 
defined from the wavenumber at which ai acquires a maximum value. It was found that 
Mp* is negative for all 9 and Moo, which indicates that ice pattern moves in the direction 
opposite to the water flow since ice grows faster just upstream of any protrusion and slower 
downstream (see FIG. 5 in Ref. |20| ). 



It is well known that a solid surface under a supercooled liquid film is morphological 
unstable, resulting in dendritic growth. The effect of the water flow on the isotherms on 
such a microscopic length scale is negligible, and the fundamental building block of the 
morphological instability of a solidification front is the Mullins-Sekerka theory. In this 
case, the amplification rate is given by al''^ = kl^,{l — {do/hQ){lth/ho)kf^{l + -ft'f)}, and the 
characteristic wavelength is Amicro = 27r{3/thC?o(l + Ki )}^^"^, where /th = i^i/V and cLq = 
TsiTCpi/L"^ are a macroscopic and microscopic characteristic length, respectively, F is the 
ice-water interface tension and Cpi is the specific heat at constant pressure of the water.— 

— 1/2 

It should be noted that Amicro depends on only Moo because /th ~ 5o ~ ^^oo , while A based 
on the macro-scale morphological instability under a supercooled liquid film herein depends 
on 9, Uoo and Q/lw, as shown in figures |5] (c), (d), (e) and (f). 
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FIG. 6. For Q/lw = 1692 [(ml/h)/cm] and x = 0.1 m, variation of dimensionless amplification rate 



cimax with (a) 9 at Uoo = 16 (km/h) and (b) 



roughness height in the experiments of Ref. 
and (d) Uoo at 6^ = 8°, for Q/ly, = 500, 1000, 1692 [(ml/h)/cm]. 
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C. Variation of amplification rate aimax with 6, u^o and Q/l 



The symbol ■ in figures |6] (a) and (b) shows the variation of measured roughness height 
by Ref. lo| with 6 for a light wind of u^o = 16 km/h and with Uoo for a mild slope of 
6 = 8°, respectively, for Q/lw = 1692 [(ml/h)/cm]. These experimental results indicate that 
roughness height increases with increasing slope and decreases with increasing wind speed. 
The variation of almax with 6 and Uoo are shown by the solid and dashed curves with the 
same parameters as used in the experiments. Small disturbances of the ice-water interface 
are assumed to be sinusoidal, expressed as 



y* = C* = 5blm[exp(a*t* + ik^x^)] = 5b{t^) sin[/c/*(a;* - VpJ^)]. 



(52) 
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In linear stability analysis, the amplitude of the ice-water interface disturbance of the most 
unstable mode increases with time as follows: 6b{t^.) = exp((Tlmaxt*)5b, where 6b = Cfc/^o 
is a dimensionless initial infinitesimal amplitude. However, the linear theory is unable to 
clarify further features related to the development of disturbances. Hence, only comparison 
of trends between the solid, dashed curves and the experimental results (■) in figures E] (a) 
and (b) is meaningful since 6b{t^) do not provide an actual roughness height of the ice- water 
interface. 

The solid curves consider the effect of the tangential and normal air shear stress distur- 
bances on the water-air interface, which are represented by and H^ in (HOj) and (HTj) . 
respectively. On the other hand, the dashed curves are obtained by neglecting and H^ in 
and (157|) . It was found that the roughness height to be expected from almax increases 
with 6 and decreases with Uoo, as shown by the solid curves in figures |6] (a) and (b), which 
show the same trends as for the experimental results (■). If the air shear stress disturbances 

(r) 

are neglected, cximax is slightly overestimated for 6, and drastically so for Mqo, as shown by 
the dashed curves in figures |6] (a) and (b), respectively. The model which takes into account 
the effect of the air shear stress disturbances on the water-air interface is supported by the 
experimental results. Streitz and Ettema state that "it is probable that roughness height 
would attain a maximum with slope, and become smaller with steeper slopes" .— However, 
for much larger values of 6, which is not shown in figure |6] (a), almax increases only grad- 
ually with 6 and never attain a maximum value. Since wind flow increases heat loss to 
air, it should be expected that the roughness height increases with increasing wind speed. 
However, the roughness height shown in figure [6] (b) is contrary to this expectation. The 
physical explanation of that will be given in the next section. Finally, figures M (c) and (d) 
show the variation of aimax with 6 at Uoo = 16 km/h and with Uoo at 6* = 8°, respectively, for 
Q/lw = 500, 1000, 1692 [(ml/h)/cm]. It was found that cr^max increases with Q/lw for any 6, 
whereas crlmax increases with Q/l^ for low Mqo, but decreases with Q/l^ for higher Uoo- 

D. Variation of the disturbed part of convective heat transfer rate at the 
water-interface with 9 and Uoo 

Since the dimensionless amplitude of the water-air interface, St = C,k/ho, and that of the 
ice-water interface, 6b = Cfc/^o are related as 6t = —fi\y^=iSb, the disturbance of the water-air 
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interface due to change in ice shape can be expressed as follows: 



y* = ^* = I + Im[(5fexp((T*t* + iki^X:^)] 

= 1 + 6b(t^)\fi\y^=i\sm[ki^{x^ - VpJ^) - 



(53) 



where \fi\y,=i y,=iY + (// |y.=i)Y^^ is the amplitude, and 6^. is the phase dif- 

(r) (i) 

ference between the water-air and ice-water interfaces. Here /; and are the real and 
imaginary parts of //. The water-air interface disturbance causes disturbance of convec- 
tive heat transfer rate from the water-air interface to the air. Its dimensionless form can 
be written as q'^^ = lm[-KadT^/dy\y=^]/ KaGa = lm[G'^^k/hoexp{aJ^ + ikux^)], where 
= Ha{ri)GaCkG^'p[<^t + ikx] is used and G'^ = (ho/So){—dHa/dri)\rf=o represents the dis- 
turbed part of the air temperature gradient at the water-air interface. It should be noted 
that q'^^ includes G'^ and ^fc. Using = —fi\y^=iCk, Qa* can be expressed as 



where \q'J = [{G'P ftX,=i - ^l..=i)^ + {G'P f[X=, + G'^ytX,=in^' is the am- 



plitude and Gg^^ is the phase difference between the disturbed heat flux at the water- 
air interface and ice-water interface. Here gT'' = {ho/6o){—dHa^ /dri)\j^=o and = 
(ho/ 6o){—dHa^ / dri)\rf=o represent the real and imaginary parts of G'^. Deflning the undis- 
turbed part of local convective heat transfer coefficient at the water-air interface and the dis- 



turbed part of it as = -KadTa/ dy\y^-hj {Tia-T^) and h'^ = lm[-KadT^/dy\y=-hJ {Tia - 



Too)], respectively, h'^/hx is equivalent to q'a^,-— 

For Q/lw = 1692 [(ml/h)/cm], figures [7] (a) and (b) show the variation of |g^^| with 6 
at Uoo = 16 km/h and with u^o at 6 = 8°, respectively. Here is evaluated from the 
wavenumber at which ai^^ acquires a maximum value. In the results represented by the solid 
curves, the effect of the tangential and normal air shear stress disturbances on the water-air 
interface is considered, which is not the case for the dashed curves. In response to the tem- 
perature distribution in the neighborhood of the ice-water interface, the amplification rate 
or roughness height of the ice-water interface is determined. Since the boundary condition 
(13^ can be written as dHi/ dy^\y^=i — G'^fi\y^=i = 0, the disturbed temperature distribution. 
Hi, in the water film is affected by q'^^. Therefore, figures E] (a) and [7] (a), as well as figures E] 
(b) and[7l^b) show the same trends with respect to 6 and Mqo, respectively. 

Figures[7](c), (d) and figures[7](e), (f) show the variation of the amplitude of the disturbed 
temperature gradient at the water-air interface, \G'J = [{Ga^^^ + {G'a^)'^Y^'^, and the ratio 



(54) 
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FIG. 7. For Q/lw = 1692 [(ml/h)/cm] and x = 0.1 m, variation of amplitude of dimensionless 
disturbed heat flux at the water-air interface, with (a) 9 at Mqo = 16 (km/h) and (b) u^o 

at = 8°. Variation of magnitude of disturbed temperature gradient at the water-air interface, 
and the ratio of amplitude of the water-air interface to that of the ice- water interface, l^^fc/Cfcl) 
with ka* for = 0° and 20° ((c) and (d)) and for Uoo = 5 and 30 (m/s) ((e) and (f)). The solid 
curves consider the effect of the tangential and normal air shear stress disturbances on the water-air 
interface, and the dashed curves do not consider this effect. 

of amplitude of the water-air interface to that of the ice-water interface, \^kKk\) against 
ka* for 9 = 0°, 20° and for Uoo = 5, 30 (m/s), respectively. For extremely small values of 
Moo, the influence of airflow on the air temperature distribution can be neglected, and 
can be approximated as (PHa/dyf = k'^^Ha- Its solution is Ha = e"'^'"*^ with the boundary 
conditions ( HT^ . which yields |G"„| = (/io/5o)^a*- For example, for Uoo = 0.01 m/s and 
Q/lw = 1692 [(ml/h)/cm], the following values are obtained: 6o = 16.1 mm, and Jiq = 0.42 
mm for 6 = 20°, Jiq = 0.56 mm for 6 = 8°. Therefore, each \G'J is represented by 0.026fca^, 
and O.OSS/cq*, which are very much below the solid curves in figures El^c) and (e). This 
indicates that the disturbed temperature gradient at the water-air interface in the presence 
of airflow becomes extremely large. 

In figure [Tl^c), IG^I for 6 = 0° is greater than that for 6 = 20°. In figure [Tl^e) , for 
Moo = 30 m/s is greater than that for Moo = 5 m/s. Therefore, it might be expected that 
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the heat transfer rate at the water-air interface for ^ = 0° and Uoo = 30 m/s is larger than 
that for 6 = 20° and Uoo = 5 m/s. However, it should be stressed that the convective heat 
transfer rate q'^^ depends not only on the disturbed temperature gradient at the water-air 
interface, G'^, but also on the amplitude of the water-air interface, = ~fi\yt=iCk- Figure [7] 
(d) shows that |^A;/Cfc| ^it 6 = 0° decreases more rapidly with ka^, than that at ^ = 20°. For 
Q/lw = 1692 [(ml/h)/cm] and Uoo = 16 km/h, ai^^ acquires a maximum value at ka* = 0.05 
for 6^ = 0° and fca^, = 0.41 for 9 = 20°. Estimating \G'J and |^fc/Cfc| at these ka*, \q'^J at 
6 = 20° is greater than at ^ = 0°, as shown in figure [7] (a). Hence, the roughness height in 
figure M (a) increases with 9. On the other hand, figure [7] (f) shows that \C,k/Ck\ decreases 
with Uoo- For Q/lw = 1692 [(ml/h)/cm] and 6 = 8°, ai^^ acquires a maximum value at 
ka* = 0.27 for Uoo = 5 m/s and ka* = 0.13 for Uoo = 30 m/s. Estimating \G'J and |^fc/Cfc| 
at these ka*, Ig'^J at Uoo = 30 m/s is less than at u^o = 5 m/s, as shown in figure [7] (b). 
Hence, the roughness height in figure [6] (b) decreases with Uoo- The undisturbed part of 
the local convective heat transfer coefficient can be written as = 0. 292 K a ^/u^^JijA^) for 
larger Uoor^ which indicates that as wind speed increases, the undisturbed part of heat 
transfer from the water-air interface to the air increases and the undisturbed ice growth 
rate, V = —h^Too/L is enhanced by the airflow. However, the disturbed part g^^ = h'^/hx 
does not necessarily increase with -Uoo- 



When the air shear stress disturbances are considered, the wavelength is A = 1.42 cm 
from fca* = 0.13 for Uqq = 30 m/s, Q/lw = 1692 [(ml/h)/cm] and 9 = 8°. If we neglect 

(r) 

the air shear stress disturbances, a* acquires a maximum value at ka* = 0.18 for the same 
parameters, and |^fc/Cfc| is overestimated for small region of fca*, as shown by the dashed 
curve in figure [7| (f). This leads to an overestimation of l^a^,!, as shown by the dashed curve 
in figure[7](b), and A = 1.03 cm from ka^. = 0.18. These results indicate that if we neglect the 
effect of the air shear stress disturbances on the water-air interface, the roughness spacing 
and height are erroneously estimated. In particular, the results show that the variation of 
the roughness height with u^o cannot predicted even qualitatively. 
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E. The most dominant term contributing to the behavior of wavelength A and 
amphfication rate al^^ 

Since the governing equations ([23]), ([24]), ([32]), ([33]), ([35]),® include 
many values: Rea, Rei, Pei, duu/ dy^\y,^=Q, duu/ dy^\y,^=i, d'^uu/dy\^=i, cos6'/Fr^, Wekf^, 
''\ ni'^ and ni'\ it is necessary to extract the most essential ones contributing to the 
behavior of A in figures [5] (c), (d) and crlmax in figures [6] (a), (b). Here, Sa \ si*^ and Ili'^'', 
Ha ■* are the real and imaginary parts of tangential and normal air shear stress disturbances, 
Tia and Ha, respectively. In the following, the values Wekf^, '^a \ '^a\ ni'^-' and na-* are 

(r) 

estimated from the wavenumber at which ai acquires a maximum value. 

First, let us consider the dependence of these values on 9. Figure [S] (a) shows that the 
variation of Rca, Rei, Pei with 6 is small except for small 6 values. For Uoo = 0.01 m/s in 
figure [8] (b), gravity-driven fiow is dominant and so, the velocity profile in the water film is 
ui* = —yl + 2?/*, which yields values dui^/ dyJ\y^=Q = 2, duu/dy^\y^=i = 0, d'^uu/dyl\y,=i = 
—2 for any 6. On the other hand, for u^o = 16 km/h in figure [8] (c), shear stress-driven 
fiow is dominant at ^ = 0° and so, the velocity profile in the water film is ui^, = which 
yields values duu/ dy^\y^=Q = 1, dui^:/ dy^\y^=i = 1, d'^ui^/ dylly^^i = 0. As 6' increases, the 
profile changes from ui* = y^ to = —yl + 2?/*. Therefore, duu/dy*\y,=o, duu/ dyj\y^=i, 
d'^uul dy1\y^=i approach the values (2, 0, -2) for larger 6 values in figure [8] (c). The effect 
of gravity and surface tension on the water-air interface is due to the terms cosO/Fr"^ = 
gho cos 9 /uf^ and Wekf^ = 'y / {piufJio)kf^ in (I37)) . respectively. The term cosO/Fr"^ for 
Uoo = 16 km/h in figure [8] (c) is greater than for u^o = 0.01 m/s in figure [8] (b), for small 
6 values. However, the difference decreases with 6, and finally the term Wekf^ is more 
dominant than the term cos 6 /Fr"^ for higher wavenumber. Figure [8] (d) shows that the 
values Tia \ \Ha ''\ and Ua^ for Uoo = 16 km/h are large for small 6 values. From (1^ 
and (I^Tl) . the values Ha \ and Ha'' decrease with 6 because Jiq decreases and uia 

increases with 6, as shown in figures HI (a) and (d). It is found from the above estimates 
that the most dominant term in the third term of 03 7p for small 6 values is cosO/Fr"^. A 
slight difference between wavelength A for Uao = 16 km/h and that for Uao = 0.01 m/s in 
figure [5] (c) appears only for very small values of 9, but the difference between them decreases 
with 9. This trend is almost the same as the variation of cos 9 /Fr"^ with 9. Since the values 
dui^/dy^\y^=i, cos 9 /Fr'^, We and IhI*^^! decrease with 9, the third term in (1371) becomes more 
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FIG. 8. For Q/lw = 1692 [(ml/h)/cm] and x = 0.1 m, variation of values (a) i?e„, Rei, Pef, (b), 
(c) duu/dy^\y,=o, duu/dy^\y^=i, d'^uu/dyl\y,=i, cos6/Fr'^, Wekf^; (d) , ni*\ with 

6. (e), (f), (g) and (h) represent the variation of these parameters with Uoo- 

effective for higher wavenumber (or equivalently, small wavelength) with 6, which results in 
a decrease of A with 6, as shown in figure [5] (c). Moreover, the amplitude of the water-air 
interface becomes smaller by the restoring force mainly due to gravity for small 6 values, 
as shown by the solid curve of 6 = 0° in figure [7| (d). As 6 increases, the action of gravity 
on the water-air interface becomes smaller since the value cos 6 /Fr'^ decreases. Hence, the 
amplitude of the water-air interface becomes larger, as shown by the solid curve of = 20° 
in figure [7] (d). This large disturbance of the water-air interface causes the increase of almax 
or roughness height with 6 in figure [6] (a). 

Second, let us consider the dependence of the above values on Uoo- Figure [8] (e) shows 
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that the value of Rca increases rapidly with Mqo, while the values of Rei and Pei increase 
only gradually. For 6 = 1° in figure [8] (f), the gravity- driven flow is dominant for small Uoo 
values, whereas the shear stress-driven flow is dominant for large Uoo values. Hence, the 
values of dui^/dy^\y^=o, dui^/ dy^\y^=i, d'^ui^/dyl\y^=i in figure E] (f) start at (2, 0, -2) and 
approach (1, 1, 0). That trend is also the case for 6' = 8° in figure E] (g) , but the shear-driven 
effect appears later than for 6 = 1°. The values cos 6 /Fr"^ and Wekf^ in figures |8] (f) and 
(g) decrease with Uoo because Jiq decreases and uia increases as Uoo increases, as shown in 
figures m (b) and (e). On the other hand, the values Tia \ Ila ■* and Ha^ in figure [8] (h) 
increase with Mqo because f HU]) and fHT]) increase with Uoo- For 6 = 1°, the term cosO/Fr'^ is 
the most dominant in the third term of fl37p for small Uoo values, but cos 6 /Fr"^ is comparable 
to Wekf^ with Uoo- On the other hand, for 6 = 8°, cosO/Fr"^ is comparable to Wekf^ for any 
Moo, as shown in figure |8] (g). The value cos 6 /Fr"^ for = 1° is much greater than that for 
^ = 8°, for small Mqo values. Hence, the variation of A with Uoo for = 1° is large compared 
to that for 9 = 8°, as shown in figure |5] (d). For Uoo = 5 m/s in figures M (g) and (h), since 
the values cos 6 /Fr'^, Wekf^, and Hq are comparable and small, the amplitude of the 
water-air interface becomes large, as shown by the solid curve of Uoo = 5 m/s in figure [7] 
(f). On the other hand, for u^o = 30 m/s in figures |8] (g) and (h), the values and H^ 
are much greater than the values cosO/Fr"^ and Wekf^. Hence, the effect of the shear stress 
disturbances on the water-air interface is more dominant than that of gravity and surface 
tension. Moreover, the amplitude of the water-air interface for small ka^, values becomes 
small compared to that for Uoo = 5 m/s, as shown by the solid curve of Uoo = 30 m/s in 
figure [7] (f). This causes the decrease of almax or roughness height with Uoo in figure [6] (b). 
As the wavenumber increases, the term Wekf^ becomes the most dominant one, and the 
water-air interface tends to be fiat because of the restoring force due to the surface tension, 
as shown in figures [7] (d) and (f). 



IV. CONCLUSION 



A theoretical model to explain the roughness characteristics in an initial aufeis (icings) 



formation observed in the experiments of Ref. 110| was proposed, from a new morpholog- 
ical instability of an ice surface during growth under a supercooled water film driven by 
gravity and air drag. A numerical method to solve complex air-water-ice multi-phase sys- 
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tern, where air and water flows and temperature flelds are highly coupled with water fllm 
thickness, was also proposed. Using linear stability analysis, roughness characteristics such 
as roughness spacing and height of the ice-water interface were derived for various water 
supply rates, plane slopes and airspeeds. Major flndings are as follows: (1) The roughness 
spacing decreases with increasing slope and airspeed, whose trends are in qualitative agree- 



ment with the experimental results of Ref. |lO|. Moreover, the roughness spacing was found 



to increase with the water supply rate. In particular, roughness spacing depends mainly on 
water layer thickness. (2) The upstream propagation of the ice-water interface disturbance 
was predicted. (3) The ampliflcation rate of the ice-water interface disturbance increases 
with slope, but decreases with airspeed. In the linear stability analysis, the ampliflcation 
rate of the ice-water interface disturbance is expected to be relevant with the initial rough- 
ness height. In the experiments,— the roughness height increased with slope, but decreased 
with airspeed. Therefore, the theoretical results herein are consistent with the experimental 
ones. Also, the ampliflcation rate was found to increase with the water supply rate for low 
airspeeds, which suggests that roughness height increases with water supply rate. The most 
important flnding of this study is that in order to predict the roughness height at higher 
airspeeds, the influence of air shear stress disturbances on the water-air interface must be 
taken into account. Otherwise, the disturbed part of the convective heat transfer rate at 
the water-air interface is not correctly predicted, and the roughness height is erroneously 
estimated. 

The roughness characteristics shown in the initial aufeis formation have common features 
with the roughness features observed in glaze icing and geological pattern formations. In 
order to extend the present model to practical aircraft and structural icing problems, we 
have to consider air, water and ice interactions near the stagnation point of objects, as 
shown in flgure[T] (c). In that case, water fllm on the object is formed by unfrozen impinging 
water droplets, the free stream velocity is not constant around the object and the angle 
changes locally along the object. Therefore, the water supply rate Q/lw used in the current 
model must be replaced by liquid water content (LWC), and the values LWC, u^o and 6 
must change locally along the position of the object. Furthermore, gravity impedes the 
supercooled water flow due to air shear stress on the upper side of the object, while it 
magnifles the effect of air shear stress on the lower side of the object. In addition to this 
asymmetry, it will be necessary to consider heat conduction into the object beneath the ice 
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sheet. Although the method and fundamental ideas developed here can be apphed to other 
phenomena by extending the current model, further research and laboratory experiments 
are necessary to validate our model proposed in this paper. 
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